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Abstract 

A study of QCD coherence is presented based on a sample of about 397 000 e + e” hadronic an¬ 
nihilation events collected at yfs = 91 GeV with the OPAL detector at LEP. The study is based 
on four recently proposed observables that are sensitive to coherence effects in the perturbative 
regime. The measurement of these observables is presented, along with a comparison with the 
predictions of different parton shower models. The models include both conventional parton 
shower models and dipole antenna models. Different ordering variables are used to investigate 
their influence on the predictions. 


1 Introduction 

Processes involving the strong interaction, described in the standard model (SM) by quantum chro¬ 
modynamics (QCD), dominate in high energy particle collisions. It is therefore important to account 
for QCD effects and to model them accurately. Colour coherence, the destructive interference effect 
between colour-connected partons, is an important aspect of high energy collisions and QCD parton 
cascades. Coherence is itself a subject of considerable interest, and QCD offers a situation in which 
coherence effects in a perturbative framework can be studied in a uniquely precise way. Furthermore, 
by testing different theoretical schemes for coherence, QCD Monte Carlo (MC) event generators (see 
Refs. [2-5] for recent reviews) can be modified to better describe the results of experiments. For 
example, in new-physics searches at the CERN LHC, QCD multijet events often represent the most 
difficult SM background to characterize. Improvements in the reliability of QCD event generators 
may help to better constrain this background. 
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The e + e _ annihilation process offers a favorable environment to study colour coherence, because 
the lack of strong interactions in the initial state allows simple and conclusive comparisons between 
experiment and theory. Previous studies of coherence in e + e _ annihilation events are presented, for 
example, in Refs. [6,7]. Within the context of a QCD shower, coherence implies an ordering condition, 
such as a requirement that each subsequent emission angle in the shower be smaller than the previous 
angle [8,9]. However, there are many ambiguities in the definition of the ordering variable and in its 
implementation. In this study, we present the first experimental tests of recently proposed [10] observ¬ 
ables designed to discriminate between coherence schemes. The data were collected with the OPAL 
detector at the CERN LEP collider at a centre-of-mass energy of y/s = 91 GeV. The observables 
examined here arc based on four-jet e + e annihilation configurations in which a soft gluon is emitted 
in the context of a three-jet topology, with two of the three jets approximately collinear. This event 
configuration has been shown to be favorable for the manifestation of coherence [11] and sensitive to 
the choice of the ordering variable in the shower [12]. 

We examine six different models for coherence, which are implemented in currently available 
QCD MC event generator programs. Specifically, we compare the default q 2 parton shower of Her- 
WIG++ [13] with angular-ordering, the p 2 j_ dip - and q^ ip -ordered dipole showers of HERWIG++ [11], 
the default p^ cvol -ordered shower of Pythia 8 [14], and the p 2 _. xn{ - and m^ nt -ordered showers of 
Vincia [15], a plugin to the Pythia 8 event generator that replaces the Pythia 8 shower with a 
shower model based on antenna functions. The definitions of the ordering variables are given in the 
next section, with more details presented in Ref. [10]. 

This paper is organized as follows. In Section 2, we define the observables to be used in the 
analysis. Sections 3, 4 and 5 present the detector, the data sample and simulation, and the data analysis, 
respectively. The results are presented in Section 6 and our conclusions in Section 7. 

2 Theoretical concepts 

2.1 Observables 

We consider hadronic events from e + e' annihilation at the Z boson peak and use the Durham kr 
clustering algorithm [16] to cluster all particles of an event into jets, keeping track of the clustering 
scales along the way. The algorithm begins by assigning all particles in an event to a list. Each entry 
in the list is called a jet. The algorithm then computes, for all pairs of four-momenta i and j in the 
event, the distance measure 


Vij = 2min (E?, E?) (1 - cos )/s , (1) 

where E t and Ej are the corresponding energies and 0,.j is the angle between objects i and j. The 
center-of-mass energy-squared is denoted by s. The pair of objects with the smallest y %! is combined 
by summing their four-momenta and the sum is added to the list while the original four-momenta are 
removed. This procedure is iterated until one entry is left in the list. To obtain an inclusive four- 
jet event sample in the perturbative regime we impose an explicit requirement on the value of the 
clustering scale at which the event goes from having four to having three jets. Denoting this scale 
(given by the value of min(y^) evaluated at the stage when the event has been clustered to four jets) 
by y 4-5-3 , we require y 4-5-3 > 0.0045 (corresponding to ln(y 4_>3 ) > —5.4), as in Ref. [10]. This 
value originates from a compromise; on the one hand a smaller value results in a data sample with 
greater statistical precision, and on the other hand a larger value provides a more direct representation 
of the shower properties. 
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Figure 1: a) The event topologies resulting from four-jet events with requirements on the angles 
between the jets: 0\2 > 2 - 7 r/ 3 , #13 > 2n/3, and 623 < 7 t/ 6 . b) Illustration of the observable 614 , the 
angle between the first and fourth jet in the latter events, c) The event topologies resulting from four- 
jet events with requirements 6 12 > 2ir/3, 0 13 > 2-k/3, 623 < 7 r/ 6 , and 624 < 7 t/ 2 . d) Illustration of 
the observable 6 * = 624 — 623 , the difference in opening angles, e) The sketch shows event topologies 
where the third and fourth clusterings occur within the same jet and hence the mass ratio M^/M^ is 
small, f) Events with large mass ratios, where the third and fourth clusterings occur on opposite sides 
of the event. 


We investigate four different observables, where for the first three we consider the event clustered 
into four jets, and order the jets in energy. To be sensitive to coherence, the angles between the jets are 
constrained such that the first (hardest) jet lies back-to-back to a nearly collinear jet pair - , formed by 
the second and third jet: 612 > 2tt/3, 6 13 > 2n/3, and 623 < vr/ 6 . The event topology resulting from 
these requirements is shown in Fig. 1 a). To investigate QCD colour coherence effects we examine 
the following observables: 

• 0 i 4 , proposed in Ref. [12]: 

The emission angle of the soft fourth jet with respect to the first jet; a sketch can be found in 
Fig. 1 b). 

• 9 *, proposed in Ref. [11]: 

A restriction on the angle between the second and fourth jet, 6 F 4 < 7 r/ 2 , is imposed in order 
to require the fourth jet to be close in angle to the nearly collinear (23) jet pan - , see Fig. 1 c). 
The observable is the difference in opening angles, 9* = 024 — 023, and is sensitive to coherent 
emission from the (23) jet system. A sketch of this observable is shown in Fig. 1 d). 

• proposed in Ref. [17]: 
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In general we have the freedom to chose the exponent ft of the 2-point energy correlation double 
ratio off ' 1 . With the choice ft = 1/5 and a nearly collinear (23) jet pair we can express the 2- 
point double ratio as C ~ E 40 \^ > d! 2 ^ 4 E v \ s /(E \with the total visible energy E v j s 
in the event. 023 4 denotes the angle between the softest jet and the (23) jet pair and analogously 
for (ft 23. Therefore the 2-point double ratio is mainly sensitive to the relative energy of the 
fourth jet. 


Strong ordering in the parton shower refers to strong ordering of the clustering scales, y 3 ^ 2 » 
y 4 —>3 w ith y( n+1 ) _m the value of the jet distance parameter in the Durham algorithm for which 

the configuration changes from n + 1 to n jets. In contrast, events with, e.g., y 4- *' 3 ~ y 3- *' 2 are more 
sensitive to the ordering condition and to situations where the same parton participates in two splitting 
processes, hence to effective 1 —> 3 splittings. For the final observable, we cluster events into two 
jets and apply the restriction y 4_>3 > 0.5 y 3 ‘ >2 . This forces events into a compressed hierarchy, i.e., a 
hierarchy without strong ordering. The investigated observable is: 

• p = proposed in Ref. [12]: 

The ratio of the invariant masses-squared of the jets at the end of the clustering process, ordered 
such that M 2 < Mj f . For “same-side” events, where one 1 —>■ 3 splitting occurs, the mass ratio 
is close or equal to zero, whereas for “opposite-side” events with 1 —> 2 <g> 1 —> 2 splittings, 
the mass ratio is larger. In Fig. 1 e) and f) we illustrate examples of these event topologies. For 
references to heavy jet masses see, e.g. Refs. [18,19]. 


To exhibit the differences between the theory models more clearly, we introduce the asymmetry 
for a given observable x, 

E n i 

N\ e f t __ i with x(i)<XQ 

bright E n i ’ 

i with x(i)>x 0 


where n, is the number of events in histogram bin i and x{i) is the bin center. The dividing point xy 
separates the regions with small and large values of x. We use this asymmetry for three of the four 


observables, 9*, EEE and p, and thus introduce three dividing points: 


Q* X~<(l/5) 
u 2 ,o ’ 


and p o- 


As in Ref. [10], we divide the full 0\4 range into three regions labelled “towards” (small $ 14 ), 
“central” (intermediate $ 14 ), and “away” (large 61 ft), denoted “T”, “C”, and “A” respectively. In the 
towards region, the first and fourth jets are collinear, while they are back-to-back in the away region. 
Events in which the fourth jet represents a wide-angle emission from the three-jet system populate the 
central region. We then consider the ratio between regions Rj and Ift : , 


E n i 

Rj iCRj 
Rk R'i 


(3) 


We define 9 different versions of the ratio, with different definitions of the regions, which are given in 
Table 1. 


2.2 Theory models 

For parton showers based on 1 —> 2 splittings of a parton I to daughters i and j, momentum conserva¬ 
tion requires that the virtuality of the branching parton must be compensated for by a recoil somewhere 
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Central/Towards Central/Away Towards/Away 


# Central region Towards region Away region Towards region 


1 

0.4 

< 

6 u /ir 

< 

0.6 

6*14/77 

< 

0.3 

#14/77 

> 

0.6 

#14/77 

< 

0.3 

2 

0.4 

< 

0 u/lT 

< 

0.6 

#14/77 

< 

0.2 

#14/77 

> 

0.7 

#14/77 

< 

0.3 

3 

0.4 

< 

On/ vr 

< 

0.6 

#14/77 

< 

0.4 

#14/77 

> 

0.8 

#14/77 

< 

0.3 

4 

0.45 

< 

#14/77 

< 

0.55 

#14/77 

< 

0.3 

#14/77 

> 

0.6 

#14/77 

< 

0.2 

5 

0.45 

< 

On/ 7 T 

< 

0.55 

#14/77 

< 

0.2 

#14/77 

> 

0.7 

#14/77 

< 

0.2 

6 

0.45 

< 

6*14/7r 

< 

0.55 

#14/77 

< 

0.4 

#14/77 

> 

0.8 

#14/77 

< 

0.2 

7 

0.35 

< 

6*14/71 

< 

0.65 

#14/77 

< 

0.3 

#14/77 

> 

0.6 

#14/77 

< 

0.4 

8 

0.35 

< 

6*14/71 

< 

0.65 

#14/77 

< 

0.2 

#14/77 

> 

0.7 

#14/77 

< 

0.4 

9 

0.35 

< 

6*14/77 

< 

0.65 

#14/77 

< 

0.4 

#14/77 

> 

0.8 

#14/77 

< 

0.4 


Table 1: Definitions of the @14 intervals used for the asymmetry ratios defined in Eq. (3). We define 
9 different versions of the ratio with the labeling of the regions given in the first column. The ratio 
between the results in the central and towards regions is based on the definitions in columns two and 
three, between the central and away regions on the definitions in columns two and four, and between 
the towards and away regions on the definitions in columns four and five. Taken from Ref. [10]. 


else in the event; we refer to parton I as the “emitter” and to the parton (system) absorbing the recoil 
as the “recoiler”. 

The six different theory models for the parton shower, mentioned in Section 1, are based on 
different formalisms and radiation functions: 

• In the collinear DGLAP formalism [20-22], each parton is evolved separately and undergoes 
1 —> 2 like branchings, which we denote pj —> piPj. In order to respect QCD coherence 
properties, a specific choice for the evolution variable [13,23], or additional vetos [24], are 
applied. The momentum balancing can either include all partons of the event, which we refer 
to as global recoils, or only one recoiler parton, which we refer to as local recoil. 

• Another formalism is based on the Catani-Seymour (CS) dipole functions [25], where a single 
parton emission from a pair of partons is considered. We denote the momenta involved in this 
splitting process with pipx ~> ViVjVk■ The full splitting probability is partitioned into two 
pieces, corresponding to partons I and K, respectively, acting as the emitter with the other 
acting as the recoiler. The recoil is limited to the longitudinal dir ection of the recoiler parton in 
the rest frame of I and K. If the dipole shower uses an evolution with ordering in transverse 
momentum, the shower correctly reproduces the soft properties of QCD. 

• In the QCD antenna (also called Lund dipoles) [26,27] picture, there is no fundamental dis¬ 
tinction between the emitter and the recoiler. Each colour-connected parton pair of an event is 
represented by an antenna and undergoes a splitting process of the form pipx —>• PiPjPk with 
a 2 —y 3 recoil prescription. A single antenna thereby accounts for the equivalent of two CS 
dipoles. 

The theory models we investigate here span all the above formalisms. For the ordering variables we 
use the notation Q] = (pt + pj) 2 , Q\ = (pj + p k ) 2 , and M] k = (p/ + p K ) 2 = (p ?; + pj + p k ) 2 , 
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for the splitting processes as stated above. For the DGLAP-based models, the parton I\ acts as the 
recoiler and can either represent a single parton (Pythia 8) or multiple partons (HERWIG++). 

In the following we briefly describe the main differences between the theory models used in this 
paper, mostly concentrated on the aspects described above. HERWIG++ q 2 [13], a parton shower 
model based on DGLAP splitting kernels, uses global recoils. The evolution is ordered in a variable 
proportional to energy times angle, 


__ 

Q 2 A M 'h< -Q 2 - Ql )' 


(4) 


The shower includes a matrix-element correction for the first emission and uses two-loop running of 
a s . The QCD coherence properties are respected due to the angular ordering of the parton branching 
cascade. The second shower model in the HERWIG++ event generator is HERWIG++ p 2 j_ dip [11], 
which is based on partitioned CS dipoles with local recoils within dipoles. The ordering variable is 
the relative transverse momentum of the splitting pair. 


2 _ Q 2 Qk( m ik ~ Q 2 ~ Qk ) 

P_Ldip — 


{M 2 - Q 2 ) 2 


(5) 


We do not apply matching or matrix-element corrections and use one-loop running of a s . The dipole 
shower with ordering in transverse momentum respects QCD coherence. As an alternative we use the 
same shower model, but with a different ordering variable. HERWIG++ q^ ip [11] orders the shower 
cascade in virtuality of the splitting pair, 


Qdip — Q 2 ■ 


( 6 ) 


As before we do not apply matching or matrix-element corrections and use one-loop running of a s . 
VINCI A p\ &nt [15] is a shower model based on antenna functions with local recoils within antennae. 
The ordering variable is the transverse momentum of the antenna, 


2 QlQl 

P± ant = -J72- 
m IK 


(7) 


Matrix-element corrections at LO [28] and NLO [29] are switched off and we use one-loop running of 
a s . Colour coherence is respected, since it is an intrinsic property of the antenna functions. Transverse 
momentum as the evolution variable is the preferred choice in VlNClA, as has been shown in Ref. [29]. 
However, we also use VlNClA m 2 nt [15] as an alternative to the transverse momentum ordering, which 
orders the shower evolution in antenna mass, defined as 


ml nt = mm(Q 2 i,Q 2 K ) . (8) 

The last shower model is Pythia 8 p 2 j_ cvo i [14], a parton shower based on DGLAP splitting kernels 
and ordered in transverse momentum, defined as 

2 _ Q 2 (Mj K - Qk)(Q 2 i + Qk) /q. 

P±evo1 (M 2 k + Q 2 ) 2 ' {) 

In contrast to the angular ordered HERWIG++ shower, local recoils within dipoles are applied. A 
matrix-element correction for the first emission is included and we use one-loop running of a s . To 
obtain QCD coherence properties, the shower applies angular vetoes. 
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Besides the shower models used in this paper, there are several other models: Ariadne [30], 
based on antenna functions, which is very similar to VlNClA; the CS dipole shower models of 
Weinzierl et al. [31], and Sherpa [32], which are similar to the HERWIG++ dipole shower; the de- 
ductor by Nagy and Soper [33], which is not interfaced with a hadronization model; and the virtuality- 
ordered final-state showers of Pythia [24,34], Nlljet [35] and Herwiri [36]. 

To compare the models on as equal a footing as possible, the shower and hadronization parameters 
have been readjusted with the PROFESSOR [37] tuning system, utilizing LEP data available through 
Rivet [38]. The tuning procedure is described in Ref. [10], where the resulting parameter values can 
also be found. 


3 OPAL experiment 

The OPAL experiment at LEP operated between August 1989 and November 2000. The detector 
components were arranged around the beam pipe, in a layered structure. A detailed description can 
be found in Refs. [39-41]. The tracking system consisted of a silicon microvertex detector, an inner 
vertex chamber, a jet chamber, and chambers outside the jet chambers to improve the precision of the 
z-coordinate 1 measurement. The jet chamber was approximately 4 m long and had an outer radius 
of about 1.85 m. This device had 24 sectors each containing 159 sense wires spaced by 1 cm. All 
tracking systems were located inside a solenoidal magnet, which provided a uniform axial magnetic 
field of 0.435 T along the beam axis. The magnet was surrounded by a lead glass electromagnetic 
calorimeter and a sampling hadron calorimeter. The electromagnetic calorimeter consisted of 11704 
lead glass blocks, divided into barrel and endcap sections, covering 98% of the solid angle. Outside 
the hadron calorimeter, the detector was surrounded by a system of muon chambers. Similar layers of 
instrumentation were located in the endcap regions. 

Since the energy resolution of the electromagnetic calorimeter is better then that of the hadron 
calorimeter, the resolution of jet directions and energies is not significantly improved by incorporat¬ 
ing hadron calorimeter information. Thus, our analysis relies exclusively on charged particle infor¬ 
mation recorded in the tracking detectors and on clusters of energy deposited in the electromagnetic 
calorimeter. 

4 Data and MC samples 

In the first phase of LEP operation, denoted LEP1 (1989 to 1995), the e + e _ center-of-mass energy was 
chosen to lie at or near the mass of the Z boson, yfs ~ 91 GeV. During the second phase of operation, 
denoted LEP2 (1995-2000), the center-of-mass energy was increased in successive steps from 130 to 
209 GeV. Interspersed at various times during the LEP2 operation, calibration runs were collected at 
the Z boson peak. In this analysis, we utilize data collected at \/s = 91.2 GeV during the the LEP2 
calibration runs. This allows us to exploit conditions when the detector was operating in its final, 
most advanced configuration. In addition, this will facilitate possible future comparisons with data 
collected under essentially identical conditions at higher energies. We use a sample corresponding to 
an integrated luminosity of 14.7 plT 1 . This sample is of sufficient size that systematic uncertainties 
dominate the statistical terms. To correct the data in order to account for experimental acceptance 
and efficiency, simulated event samples produced with MC event generators are used. The process 

1 OPAL uses the right-handed coordinate system defined with the a;-axis pointing towards the center of the LEP ring, the 
positive z points along the direction of the e~ beam and the y -axis upwards, r is the coordinate normal to the beam axis 
and the polar angle 9 and the azimuthal angle ifi are defined with respect to x and a. 
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e + e~ —> qq is simulated using Pythia 6.1 [42] at y/s = 91.2 GeV. Corresponding samples using 
Herwig 6.2 [43,44] are used for systematic checks. We examine the MC events at two levels. 
We refer to “hadron level” as events without event selection, and without simulation of the detector 
acceptance and resolution, for which all particles with lifetimes less than 300 ps decay. In contrast, 
“detector level” refers to MC events that are processed through the GEANT4-based [45] simulation 
of the OPAL detector, called Gopal [46], and that have been reconstructed using the same software 
procedures that are applied to the data. The MC events generated for the detector-level samples are the 
same as the hadron-level samples except that Kg mesons and weakly decaying hyperons are declared 
to be stable, as these particles can interact with detector material before decaying, and so their decays 
arc handled within the Geant framework. 

In addition, for comparisons with the corrected data, large samples of hadron-level MC events 
are employed, using the event generators Herwig++ 2.7.0 [47], Pythia 8.176 [48], and Vin- 
CIA 1.1.0 [15] interfaced with the hadronization model of Pythia 8.176. 

5 Data analysis 

5.1 Selection of events 

The same criteria for the selection of charged tracks and electromagnetic clusters are applied as de¬ 
scribed in Ref. [49]. Charged tracks are required to have transverse momentum relative to the beam 
axis larger than 0.15 GeV, and photons to have energies larger than 0.10 GeV (0.25 GeV) in the bar¬ 
rel (endcap) region of the electromagnetic calorimeter. The selection of hadronic annihilation events 
is the same as described in Ref. [50]. Briefly, a minimum of five charged Lacks is required, and a con¬ 
tainment condition | cos$t| < 0.90 is applied, where 9 x is the polar angle of the thrust axis [51,52] 
with respect to the beam axis, calculated using all accepted charged tracks and electromagnetic clus¬ 
ters. A total of 397 452 candidate hadronic annihilation events are selected, with a negligible expected 
background. 

Since the energy loss due to initial-state radiation is highly suppressed at the Z peak, we do not 
apply a cut to that effect. However, radiative corrections are applied by requiring yfs — \fs' < 1 GeV 
for the MC detector-level samples used to correct the data, where y/s* is the effective center-of-mass 
energy after initial-state radiation. 

5.2 Reconstruction and correction 

For each of the accepted events, the values of all observables described in Section 1 are computed. 
To avoid double-counting of energy between tracks and electromagnetic clusters, an energy-flow al¬ 
gorithm [53,54] is applied, which matches the tracks and clusters and retains only those clusters that 
are not associated with a track. 

Figure 2 shows a comparison of the uncorrected data with the detector-level predictions of Pythia 6 
and Herwig 6 for the # 14 , p, 9 *, and variables. The 9\ 4 and 9* variables are normalized by a 

factor of 7T. The simulations are seen to provide a generally adequate description of the measurements. 

To correct the data for detector and resolution effects, we implement an unfolding procedure 
based on the RooUnfold [55] framework. We use the iterative Bayes method, as proposed by 
D’Agostini [56], with four iterations, which is the recommendation from Ref. [55]. A necessary in¬ 
gredient for the unfolding is the response matrix of the MC event generator used for the correction 
procedure. For the standard analysis, Pythia 6 is used to determine the response matrix. The re¬ 
sponse matrix gives the bin-to-bin migration from the hadron to the detector level, and vice versa. 


a) Angle between 1st and 4th jet, d 14 /7r 


b) Ratio of jet masses, p 




Figure 2: The uncorrected distributions of a) the emission angle 6 14 , b) the mass ratio p = Alj/Mf f , 
c) the difference in opening angles 9*, and cl) the 2 -point double ratio in comparison with the 

predictions of the Herwig 6 and Pythia 6 Monte Carlo event generators at the detector level. The 
error bars indicate the statistical uncertainties. 


9 













In order to obtain reliable results for the corrected distributions, we adjust the bin widths of the his¬ 
tograms such that the probability for a hadron-level event to migrate to a different bin at the detector 
level is less than 50%. 

The corrected distributions are presented in Fig. 3 and Tables 2-3. Tables 2-3 include the co- 
variance matrices calculated with RooUnfold. The statistical uncertainties are given by the square 
root of the corresponding diagonal element in the covariance matrices. Systematic uncertainties are 
discussed in Section 5.3. Figure 3 includes the predictions of Pythia 6 and Herwig 6 at the hadron 
level. The differences between the MC predictions and the data are seen to be similar to those ob¬ 
served at the detector level (Fig. 2), demonstrating that the correction procedure does not introduce a 
discernible bias. 

The values of the derived distributions, i.e., the ratios of the different regions for 0\4 and the 
asymmetry for the other observables, are listed in Table 4. The quantities are determined by summing 
and dividing the histogram entries. The statistical uncertainties are evaluated from propagation of 
errors, while the systematic uncertainties are determined as described in Section 5.3. 

5.3 Systematic uncertainties 

Systematic uncertainties are evaluated by repeating the analysis with different selection requirements 
and with variations in the correction procedure. Specifically, we consider the following: 

• The requirement on the thrust angle direction is changed to | cos#t| < 0.7 from the default 
| cos 9t\ < 0.9. 

• The minimum number of charged tracks is increased to seven from the default of five. 

• Variation of the reconstruction procedure: All tracks and clusters are taken into account. In this 
case the detector correction takes care of the double counting. 

• Herwig 6 is used in place of Pythia 6 to determine the response matrix. 

The systematic uncertainty is determined for each variation from the bin-by-bin difference in the 
corrected distributions with respect to the standard result. The total systematic uncertainty is given by 
the quadrature sum of the individual terms. The total uncertainty of the data is defined by summing 
the statistical and systematic contributions in quadrature. 

As additional systematic checks on the unfolding procedure, we consider the following variations: 

• We use the unfolding method with three and five instead of four iterations. 

• Instead of the iterative method, we use the singular value decomposition, as proposed by Hocker 
and Kartvelishvili [57] and implemented in RooUnfold. 

We find the systematic variations that arise from these two checks to be smaller or comparable to 
the variation observed when using Herwig 6 in place of Pythia 6. Since adding all these effects 
together would likely double count the uncertainty associated with the unfolding procedure we do not 
add the observed differences to the systematic uncertainty. 

6 Comparison with Monte Carlo models 

In this section, we present a comparison between the coherence schemes described in Section 2.2 and 
the data. For this purpose, samples of 5 x 10 6 events are generated for each MC model, using the 
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Figure 3: The corrected distributions of a) the emission angle # 14 , b) the mass ratio p = M'j/M'fj, 
c) the difference in opening angles 6 *, and cl) the 2 -point double ratio in comparison with 

the predictions of the Herwig 6 and Pythia 6 Monte Carlo event generators at the hadron level. 
The error bars limited by the horizontal lines indicate the statistical uncertainties, while the total 
uncertainties correspond to the full error bars. 
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du/'K a 1 da/d(9*/w) 


0.00-0.15 
0.15-0.20 
0.20 - 0.25 
0.25 - 0.30 
0.30 - 0.35 
0.35 - 0.40 
0.40 - 0.45 
0.45 - 0.50 


0.4631 ±0.0454 ±0.1670 
1.6474 ±0.1345 ±0.4193 
1.4728 ±0.1172 ± 0.3055 
1.5833 ±0.1139 ± 0.0664 
1.5249 ± 0.1004 ± 0.2382 
1.6383 ±0.1109 ±0.1032 
1.5172 ± 0.0964 ±0.1741 
1.6025 ± 0.1053 ± 0.0624 


0.50-0.55 
0.55 - 0.60 
0.60 - 0.65 
0.65 - 0.70 
0.70 - 0.75 
0.75 - 0.80 
0.80 - 0.85 
0.85-1.00 


1.6381 ±0.1000 ± 0.0557 
1.4319 ±0.1024 ±0.1045 
1.2758 ± 0.0945 ±0.1736 
0.9020 ± 0.0708 ± 0.0741 
0.7668 ± 0.0648 ±0.1853 
0.8999 ±0.0798 ± 0.2725 
0.5220 ±0.0616 ± 0.0566 
0.0626 ± 0.0087 ± 0.0223 


Correlation Matrix 
1.000 


0.252 

-0.037 

-0.048 

-0.068 

-0.055 

1.000 

0.105 

-0.062 

-0.081 

-0.064 

1.000 

0.103 

-0.108 

-0.087 

1.000 

0.190 

-0.082 

1.000 

0.141 

1.000 


-0.068 

-0.077 

-0.082 

-0.118 

-0.106 

0.164 

1.000 

-0.055 

-0.072 

-0.084 

-0.081 

-0.114 

-0.105 

0.153 

-0.073 

-0.071 

-0.099 

-0.093 

-0.113 

-0.104 

-0.123 

-0.068 

-0.068 

-0.078 

-0.078 

-0.088 

-0.085 

-0.111 

-0.072 

-0.065 

-0.083 

-0.086 

-0.097 

-0.089 

-0.095 

-0.066 

-0.061 

-0.092 

-0.096 

-0.091 

-0.062 

-0.096 

-0.050 

-0.074 

-0.092 

-0.088 

-0.093 

-0.038 

- 0.112 

-0.073 

-0.084 

-0.091 

-0.091 

-0.044 

-0.116 

-0.129 

-0.075 

0.003 

-0.038 

-0.073 

-0.105 

-0.078 

-0.103 

-0.072 

0.006 

-0.091 

-0.092 

-0.088 

-0.091 

-0.090 


1.000 

0.273 

- 0.100 

-0.111 

-0.094 

-0.108 

1.000 

0.292 

-0.090 

-0.130 

-0.098 

1.000 

0.181 

-0.097 

-0.109 

1.000 

0.290 

- 0.110 

1.000 

0.273 

1.000 


-0.099 

-0.107 

- 0.110 

-0.116 

-0.108 

0.259 

1.000 

-0.099 

-0.087 

-0.073 

-0.096 

-0.090 

-0.079 

0.179 

-0.045 

-0.094 

-0.031 

-0.095 

-0.104 

-0.092 

0.055 


1.000 

0.190 


1.000 


Table 2: The normalized corrected data and the correlation matrix at the hadron level for the emission angle 6 * 14 . The first uncertainty is statistical 
and the second systematic. 


p 


Correlation Matrix 


a 1 da/dp 

0.00 -0.06 1.0108 ± 0.0589 ±0.1137 

0.06 -0.15 1.1474 ±0.0470 ± 0.0934 

0.15 -0.38 0.7278 ±0.0190 ± 0.0308 

0.38 - 0.69 1.0901 ± 0.0151 ± 0.0344 

0.69 - 1.00 1.0670 ± 0.0207 ± 0.0254 


1.000 

-0.007 

-0.530 

-0.356 


1.000 

0.088 1.000 
-0.490 0.387 


9*1 

7r 

a 1 d<r/d(f?*/ 7r ) 

Correlation Matrix 





-0.05- 

-0.04 

0.2849 ±0.0497 ±0.1101 

1.000 






0.04- 

0.10 

1.8538 ± 0.1850 ± 0.2395 

-0.013 

1.000 





0.10- 

0.16 

2.1298 ± 0.1579 ±0.3549 

0.358 

0.367 

1.000 




0.16 - 

0.22 

2.7726 ±0.1791 ± 0.5958 

-0.188 

0.010 

0.232 

1.000 



0.22 - 

0.28 

3.6151 ±0.2158 ±0.2605 

-0.250 

-0.327 

-0.092 

0.278 

1.000 


0.28 - 

0.34 

4.1203 ±0.2477 ±0.2650 

-0.219 

-0.457 

-0.408 

-0.295 

0.245 

1.000 

0.34- 

0.43 

1.1652 ±0.1110 ±0.1946 

-0.191 

-0.335 

-0.376 

-0.382 

-0.168 

0.304 1.000 










c ( 1/5) 

a” 1 dtr/dCf /5) 

Correlation Matrix 





0.355 - 

0.402 

4.3798 ± 0.1680 ±0.1447 

1.000 






0.402 - 

0.424 

11.2510 ± 0.2603 ±0.7473 

0.314 

1.000 





0.424 - 

0.446 

11.0688 ± 0.2557 ±0.1536 

-0.385 

0.218 

1.000 




0.446 - 

0.468 

11.8690 ± 0.3771 ±0.7888 

-0.458 

-0.493 

0.187 

1.000 



0.468 - 

0.495 

1.5553 ±0.1137 ± 0.2591 

-0.260 

-0.387 

-0.225 

0.299 

1.000 



Table 3: The normalized corrected data and the correlation matrix at the hadron level for the mass ratio p = Mj/Alf r , the difference in emission 
angles 9*, and the 2-point double ratio C^ 5 \ The first uncertainty is statistical and the second systematic. 


# Central/Towards 


1 1.0159 ± 0.0535 ±0.1059 

2 2.0383 ±0.1447 ± 0.7756 

3 0.6687 ± 0.0304 ± 0.0271 

4 0.5319 ± 0.0324 ±0.0432 

5 1.0671 ± 0.0825 ± 0.3804 

6 0.3501 ±0.0192 ±0.0139 

7 1.4942 ±0.0740 ±0.1579 

8 2.9979 ± 0.2060 ± 1.1161 

9 0.9836 ±0.0411 ± 0.0339 


6* 

/-bright 

0.10 

0.1586 ±0.0147 ± 0.0268 

0.16 

0.3599 ±0.0240 ±0.0179 

0.22 

0.7575 ±0.0441 ±0.1196 

0.28 

1.8402 ±0.1137 ±0.2619 


# 

Central/Away 

# 

Towards/Away 

1 

1.3591 ± 0.0675 ±0.0939 

1 

1.3378 ±0.0745 ±0.2070 

2 

2.6045 ±0.1591 ±0.4134 

2 

2.5637 ±0.1695 ±0.5789 

3 

8.7205 ± 0.8699 ± 1.2979 

3 

8.5840 ±0.8834 ± 2.0061 

4 

0.7116 ±0.0415 ±0.0618 

4 

0.6668 ± 0.0489 ± 0.2315 

5 

1.3636 ± 0.0932 ± 0.2226 

5 

1.2778 ±0.1041 ±0.5149 

6 

4.5656 ± 0.4765 ± 0.7924 

6 

4.2784 ± 0.4851 ± 1.7585 

7 

1.9989 ± 0.0926 ±0.1683 

7 

2.0323 ±0.0995 ±0.2145 

8 

3.8307 ± 0.2238 ± 0.6735 

8 

3.8947 ±0.2358 ±0.7391 

9 

12.8260 ± 1.2589 ± 1.9634 

9 13.0404 ± 1.2966 ± 2.3652 


Po 

^"left/-bright 

r { 1/5) 
°2,0 

Meft /-bright 

0.15 

0.1961 ±0.0069 ±0.0212 

0.42 

0.8294 ±0.0239 ±0.0503 

0.38 

0.4955 ±0.0121 ±0.0312 

0.45 

2.2991 ±0.0767 ±0.1346 


Table 4: The corrected data for the derived distributions. The upper three tables list the results for asymmetry ratios defined in Eq. (3), with 
the definitions of the towards, central, and away regions given in Table 1. The bottom three tables list the results for the asymmetries defined for 
the other observables. The first uncertainty is statistical and the second systematic 


HERWIG++ q 2 HERWIG++ p 2 ±dip HERWIG++^ ip PYTHIA 8 p 2 ±evol VlNCIA Plant VlNCIA m 2 nt 


6»i 4 /7T 

4.8 

(9.968 

X 

10- 1 ) 

12.8 

(6.850 

X 

10" 1 ) 

16.4 

(4.230 

X 

10" 1 ) 

8.3 

(9.401 

X 

10" 1 ) 

8.1 

(9.718 

X 

10" 1 ) 

7.1 

(9.444 

X 

io- 

- 1 ) 

9 */tt 

6.1 

(5.243 

X 

10" 1 ) 

2.8 

(9.032 

X 

10" 1 ) 

8.8 

(2.683 

X 

10" 1 ) 

2.5 

(9.239 

X 

10" 1 ) 

3.3 

(6.496 

X 

10" 1 ) 

5.1 

(8.552 

X 

io- 

- 1 ) 

C 0-/V 

16.2 

(6.196 

X 

io~ 3 ) 

6.0 

(3.083 

X 

10" 1 ) 

41.2 

(8.498 

X 

10 -8 ) 

4.3 

(5.090 

X 

10" 1 ) 

5.3 

(8.013 

X 

10" 1 ) 

2.3 

(3.753 

X 

io- 

- 1 ) 

P 

6.6 

(2.496 

X 

10" 1 ) 

10.1 (7.15 

X 

10" 2 ) 

42.8 

(4.098 

X 

10" 8 ) 

2.6 

(7.633 

X 

10" 1 ) 

6.8 

(5.682 

X 

10" 1 ) 

3.9 

(2.354 

X 

io- 

- 1 ) 

( V ) 


(4.442 

X 

10” 1 ) 


(4.92 

X 

10" 1 ) 


(1.728 

X 

10" 1 ) 


(7.841 

X 

10" 1 ) 


(7.477 

X 

io- 1 ) 


(6.026 

X 

io- 

- 1 ) 


Table 5: The \' 2 values (with p-values in parentheses) for the four observables considered in this analysis. The last line gives the average p-value, 
(p), for the four observables. 























tuned parameter sets mentioned in Section 2.2. We present the predictions for the different schemes 
in terms of the observables defined in Section 1. As a measure of the level of agreement with data, we 
calculate the significance, defined as 


MCj - Vi 
= - 

VVi 


( 10 ) 


where MC, and V, represent the predicted and observed values in bin i of a distribution, with ax> t the 
corresponding uncertainty in Vi. In the following, we present a distribution of the significance in a 
plot below the distribution of the variables. In addition we calculate the ;y 2 values for the distribution 
as 


Wins 

x 2 = (MC - V) T y - l {MC -V)= ( MC - V) t {V- x ) l3 {MC - V)j , (11) 

i,j =1 


where V is the full covariance matrix, representing statistical terms, with systematic uncertainties 
added to the diagonal elements. We present the x 2 results and cones ponding p-values in Table 5. The 
latter are calculated with the ROOT [58] program and give the probability that the deviations of the 
MC predictions from the data are consistent with the evaluated uncertainties. 


6.1 Angle between first and fourth jet: 0 14 

In Figs. 4 a) and b) we show the normalized distribution of the emission angle of the soft fourth jet 
from the nearly collinear three-jet system, 6 * 14 . All models are found to provide adequate descriptions 
of the data, except that the HERWIG++ p 2 j_ d[p model lies about three standard deviations above the 
measurements for a narrow region around 0 44 ~ 0.77T. Even for this model, the p-value is 0.42 
(Table 5), implying overall compatibility with the experimental results. 

We show the ratio C/T of the central-to-towards regions, which gives the relative amount of wide- 
angle to collinear emissions, in Figs. 5 a) and b). For the p^ dip -ordered dipole shower of HERWIG++ 
and the parton shower of Pythia 8 we find nearly perfect agreement with the data for all nine C/T 
regions (Table 1). The HERWIG++ q 2 model and the VlNClA p 2 _ iUlt model lie below the data by up to 
two standard deviations in some regions, while the VlNClA m 2 nt model lies about two standard devi¬ 
ations above the data in all regions. The two Vincia models exhibit the expected behavior: When the 
antenna mass is used as the evolution variable, soft wide-angle emissions are preferred over collinear 
ones, which leads to higher values for the relative level of wide-angle to collinear emissions. This 
demonstrates the sensitivity of the 6 \4 variable to the choice of evolution scheme. The largest devia¬ 
tion from the data in Figs. 5 a) and b) is observed for the q dip -ordered dipole shower of HERWIG++, 
for which the predictions lie up to around three standard deviations below the data in some regions. 
Thus, this model predicts too many collinear emissions compared to wide-angle emissions. 

In Figs. 5 c) and d) we show a comparison of the MC predictions to the data for the ratio C/A 
of the central-to-away regions. This ratio measures the relative amount of wide-angle emissions to 
emissions in a backwards direction, away from the leading jet and near to the collinear (23) jet pair. For 
the HERWIG++ q 2 model and for VlNClA, we find a good agreement with the data and observe small 
differences for the different evolution variables of VINCI A. Pythia 8 and the p^_ dip -ordered dipole 
shower of F1ERWIG++ lie below the data, by around one and two standard deviations, respectively, and 
thus predict too few wide-angle emissions compared to the backwards emissions. The HERWIG++ 
q dip model lies around one standard deviation above the data and thus predicts relatively too many 
wide-angle emissions. The observations the ratios C/T and C/A are confirmed by the measurements 
of the ratio T/A of the towards-to-away regions, presented in Figs. 5 e) and/). 
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a) Angle between 1st and 4th jet, 0 14 , Herwig++ b) Angle between 1st and 4th jet, d 14 , Vincia, Pythia 8 




9 14 /tt 6 u /n 

Figure 4: The corrected distribution of the emission angle 9 14 of the soft fourth jet in comparison 
with the predictions of a) Herwig++ and b) Pythia 8 and Vincia. The thin solid lines correspond 
to HERWIG++ with angular-ordering (q 2 ), the thick solid lines to the dipole shower of HERWIG++ 
with ordering in p^ dip , and the dash-dotted lines to ordering in q dip . VINCIA with ordering in p^_ ant 
is shown with medium solid lines, ordering in m 2 nt with dashed lines, and PYTHIA 8 is shown with 
dotted lines. The error bars limited by the horizontal lines indicate the statistical uncertainties, while 
the total uncertainties correspond to the full error bars. The ratio plots show the deviation of the 
predictions from the data in units of the total uncertainty. 


6.2 Difference in opening angles: 9 * 


In Figs. 6 a) and b) we show the normalized distribution of the difference in opening angles between 
the third and the fourth jet with respect to the second jet, 9* = #24 — $ 23 - All models are seen to 
provide an adequate description of the data, with the exception of the region around 9* ss 0.077T 
(second bin of Figs. 6 a) and b)), where the models predict somewhat fewer events than are observed. 
The largest discrepancy in this region arises from the HERWIG++ q dip model. 

We show the asymmetry as a function of the dividing point (f {] in the Figs. 6 c) and cl). The largest 
discriminating power is found for 9q* = 0.16-7T, where the ry cP|) -ordered dipole shower of HERWIG++ 
generates a deviation of almost four standard deviations with respect to the data. The number of 
events with large differences in the opening angles of the third and fourth jets is overestimated by 
this non-coherent shower model. The jr dip -ordered HERWIG++ shower, based on the same shower 
kernels, but respecting coherence due to the choice of evolution variable, gives a better description of 
the asymmetry. This emphasizes the need for coherence in order to describe the data properly. 
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a) Central/Towards for 8 U , Herwig++ 



b) Central/Towards for 8 U , Vincia, Pythia 8 



c) Central/Away for 8 u , Herwig++ 





f) Towards/Away for 8 u , Vincia, Pythia 8 



Figure 5: The corrected data for the derived distributions in comparison with the predictions of a) 
Herwig++ and b) Pythia 8 and Vincia. The thin solid lines correspond to Herwig++ with 
angular-ordering (q 2 ). the thick solid lines to the dipole shower of HERWIG++ with ordering in 
and the dash-dotted lines to ordering in q 2 ]]p . VINCIA with ordering in p 2 _. iUi is shown with medium 
solid lines, ordering in m 2 nt with dashed lines, and Pythia 8 is shown with dotted lines. The error 
bars limited by the horizontal lines indicate the statistical uncertainties, while the total uncertainties 
correspond to the full error bars. The ratio plots show the deviation of the predictions from the data in 
units of the total uncertainty. 
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a) Difference in opening angles, 9*, Herwig++ b) Difference in opening angles, 9*, Vincia, Pythia 8 



c) Asymmetry for 9 *, Herwig++ d) Asymmetry for 9 *, Vincia, Pythia 8 



Figure 6: The distribution of the difference in opening angles 0* for a) HERWIG++ and b) Pythia 8 
and Vincia. The asymmetry with respect to the dividing point 9q is shown for c) HERWIG++ and 
d) Pythia 8 and Vincia. The thin solid lines correspond to FIerwig++ with angular-ordering (q 2 ), 
the thick solid lines to the dipole shower of HERWIG++ with ordering in p^_ dip , and the dash-dotted 
lines to ordering in q 2 ]]p . VINCIA with ordering in q> 2 _. uit is shown with medium solid lines, ordering 
in m 2 nt with dashed lines and Pythia 8 is shown with dotted lines. The error bars limited by the 
horizontal lines indicate the statistical uncertainties, while the total uncertainties correspond to the 
full error bars. The ratio plots show the deviation of the predictions from the data in units of the total 
uncertainty. 
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6.3 2-Point double ratio: C^ 1 ^ 

For the normalized distribution of the 2-point double ratio, shown in Figs, la) and b), we 

find rather large deviations between the data and the MC prediction for most of the shower models. 
We again find that the (/j ip -ordered HERWIG++ shower exhibits the largest discrepancies. Only two 
models, Pythia 8 and the VlNClA m ant model, yield p-values larger than 50% (Table 5). 

In Figs. 7 c) and d), we show the asymmetry in the C 2 ' ' ^ variable as a function of the dividing 
point C 2 ^ l] ■ Since C^ T>) is proportional to the energy E\ of the fourth jet, the asymmetry in C 2 * ^ 
measures the relative number of events of soft versus hard fourth-jet emissions. We observe large 
deviations from the data, at the level of four standard deviations, for the HERWIG++ q 2 i]p model, 
which underpredicts the relative fraction of events with a very soft fourth jet. A similar discrepancy, 
at the level of around 2.5 standard deviations, is observed for the HERWIG++ q 2 model. The two 
versions of VINCI A exhibit deviations of about one standard deviation in the opposite sense, i.e., 
VlNClA m 2 nt somewhat overpredicts the level of hard fourth-jet emissions, whereas VlNClA p]_ ant 
predicts too few hard fourth-jet emissions. In contrast, the HERWIG++ p^_ dip and Pythia 8 models 
are in nearly perfect agreement with the data. 

6.4 Mass ratio: p = Mf /Mf r 

The normalized distributions of the p = variable are shown in Figs. 8 a) and b). For 

the Pythia 8 and the two VlNClA models, we find reasonable overall agreement with the data, with 
differences on the level of two standard deviations or less. The HERWIG++ models demonstrate larger 
differences, with discrepancies reaching the level of four standard deviations for the HERWIG++ g dip 
model. 

In Figs. 8 c) and d) we show the asymmetry of the p variable as a function of the dividing point pq. 
This asymmetry is sensitive to the relative number of same-side versus opposite-side events, whose 
definitions were given in Section 2.1. This asymmetry is seen to provide discrimination between most 
of the shower models. The Pythia 8 and Herwig++ q 2 models yield predictions that lie within one 
standard deviation of the data. However, the HERWIG++ q 2 - p model predicts too small an asymmetry 
by about four standard deviations, meaning that there are too few same-side compared to opposite- 
side events. The two VlNClA models also predict too few same-side events, but only at the level of 
around one standard deviation. In contrast, the HERWIG++ p^ dip model predicts relatively too many 
same-side events, at the level of two standard deviations. 


7 Summary and conclusion 

We have presented measurements of distributions in e + e~ annihilations at y/s = 91.2 GeV that are 
sensitive to QCD colour coherence, the ordering parameter in parton showers, and to whether four-jet 
events arise from two separate 1 —I 2 splittings or from a 1 —> 3 splitting. The data, corresponding 
to a sample of about 397 000 hadronic annihilation events, were collected with the OPAL detector at 
LEP. The event selection criteria are defined in a way to minimize the influence of non-perturbative 
(hadronization) effects. We compared the data with six different models for the parton shower, based 
on the HERWIG++, Pythia 8 , and VlNClA Monte Carlo event generator programs, which differ in 
the choice of the radiation function, ordering variable, and recoil strategy. Each of the six models was 
found to be in general agreement with the data. However, interesting differences between the models 
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a) 2-point double ratio, C ^'", Herwig++ b) 2-point double ratio, Vincia, Pythia 8 



c) Asymmetry for cf^ 5) , Herwig++ 


d) Asymmetry for Vincia, Pythia 8 



Figure 7: The distribution of the difference in opening angles for a) HERWIG++ and b) 

Pythia 8 and Vincia. The asymmetry with respect to the dividing point C'^'q 5 " 1 is shown for c) 
Herwig++ and d) Pythia 8 and Vincia. The thin solid lines correspond to Herwig++ with 
angular-ordering ( q 2 ), the thick solid lines to the dipole shower of HERWIG++ with ordering in p^ dip , 
and the dash-dotted lines to ordering in g^ ip . VINCIA with ordering in p 2 _, ini is shown with medium 
solid lines, ordering in m 2 nt with dashed lines and Pythia 8 is shown with dotted lines. The error 
bars limited by the horizontal lines indicate the statistical uncertainties, while the total uncertainties 
correspond to the full error bars. The ratio plots show the deviation of the predictions from the data in 
units of the total uncertainty. 
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a) Ratio of jet masses, p, Herwig++ b) Ratio of jet masses, p, Vincia, Pythia 8 



c) Asymmetry for p, Herwig++ d) Asymmetry for p, Vincia, Pythia 8 


10 ° 


2 

0 

-2 

-4 


Po 

Figure 8: The distribution of the difference in opening angles p = Mf/M 2 r for a) HERWIG++ 
and b) Pythia 8 and Vincia. The asymmetry with respect to the dividing point po is shown for 
c) Herwig++ and d) Pythia 8 and Vincia. The thin solid lines correspond to Herwig++ with 
angular-ordering (q 2 ), the thick solid lines to the dipole shower of HERWIG++ with ordering in p5_ dip , 
and the dash-dotted lines to ordering in q 2 y ]p . VINCIA with ordering in p 2 _. ini is shown with medium 
solid lines, ordering in m 2 nt with dashed lines and Pythia 8 is shown with dotted lines. The error 
bars limited by the horizontal lines indicate the statistical uncertainties, while the total uncertainties 
correspond to the full error bars. The ratio plots show the deviation of the predictions from the data in 
units of the total uncertainty. 
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and between some of the models and the data were observed when asymmetries in the distributions 
were examined. 

Until now it was nearly impossible to distinguish between the predictions of Pythia 8 and VlN- 
CIA, or between the different variants of VlNClA. Our study of the asymmetry of the ratio of squared 
jet masses, shown in Fig. 8 d), shows that VlNClA predicts somewhat too many opposite-side events 
(i.e., events with two 1 —I 2 splittings) compared to same-side events (i.e., events with a 1 —y 3 
splitting), and that the data prefer Pythia 8. We find that the different variants of VlNClA can be 
distinguished using the central-to-towards (Fig. 5 b)) and central-to-away (Fig. 5 d)) ratios in the #14 
variable, which indicate that the VlNClA variant based on antenna mass-squared evolution predicts 
somewhat too many wide-angle emissions for the soft fourth jet, compared to collinear emissions. 

To summarize the results of our study, we show the average p-value for each model, calculated 
from the four values of the single observables, in the bottom row of Table 5 . The variant of HERWIG++ 
with a -ordered dipole shower is found to provide the least satisfactory description of the data. This 
model does not contain coherence; it has intentionally been introduced to confront it with coherent 
evolution. Thus our results emphasize the importance of incorporating coherence into the description 
of the QCD multijet process. Since Herwig++ uses the cluster [23] and Pythia 8 and Vincia use 
the Lund string [59, 60] hadronization model, a direct comparison of the predictions from the two 
groups of shower models is somewhat ambiguous. It would be interesting to perform a comparison 
based on use of the same hadronization model for all models. However, when comparing all shower 
models together, we find PYTHIA 8 and VlNClA with evolution in transverse momentum to give the 
best description of the measurements presented here. 
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